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Abstract 

The  problem  of  parameter  estimation  is  considered  for  the  case  of  mathematical  models  for  polymer  electrolyte  membrane  fuel  cells 
(PEMFCs).  An  algorithm  for  nonlinear  least  squares  constrained  by  partial  differential  equations  is  defined  and  applied  to  estimate  effective 
membrane  conductivity,  exchange  current  densities  and  oxygen  diffusion  coefficients  in  a  one-dimensional  PEMFC  model  for  transport  in 
the  principal  direction  of  current  flow.  Experimental  polarization  curves  are  fitted  for  conventional  and  low  current  density  PEMFCs.  Use  of 
adaptive  mesh  refinement  is  demonstrated  to  increase  the  computational  efficiency. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Mathematical  models  are  now  widely  used  to  simulate  the 
behavior  of  fuel  cell  devices.  Such  models  are  desirable  be¬ 
cause  they  can  provide  general  trends  as  well  as  quantitative 
measures  of  relative  changes  in  performance  for  the  device 
as  model  parameters  are  varied.  The  models  can  also  provide 
detailed,  localized  data  that  are  frequently  unavailable  from 
experimental  measurements  within  an  operating  cell. 

However,  the  detailed  data  and  parametric  studies  pro¬ 
duced  by  computer  simulations  should  not  be  accepted  with¬ 
out  appropriate  validation  with  experimental  data.  This  vali¬ 
dation  may  include  reduction  of  modeling  error  by  comparing 
a  chosen  mathematical  model  to  other  potential  models  or  re¬ 
duction  of  parameter  error  by  adjusting  the  parameters  of  the 
chosen  model.  In  each  case,  the  goal  is  to  minimize  the  dis¬ 
crepancy  between  the  computed  simulation  and  the  observ¬ 
able  experimental  data.  Estimation  of  the  model  parameters 
is  the  focus  of  this  study. 

The  problem  of  parameter  estimation  can  be  summarized 
in  general  terms  as  follows.  Suppose  that  we  have  developed  a 
mathematical  model  which  takes  the  form  of  a  system  of  par¬ 
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tial  differential  equations  (PDEs)  for  the  solution  variables, 
as  well  as  a  set  of  model  parameters,  both  discrete  (finite)  or 
spatially  distributed  (infinite).  Suppose  also  that  certain  ex¬ 
perimental  data  are  available,  which  we  assume  correspond 
to  evaluations  of  various  observable  outputs  of  the  solution 
variables.  The  problem  is  then  to  find  an  optimal  set  of  pa¬ 
rameter  values  that  minimizes  the  error  between  the  observed 
output  functionals  of  the  experimental  data  and  those  outputs 
calculated  using  the  solution  to  the  mathematical  model. 

In  the  case  of  fuel  cells,  while  there  have  been  some  re¬ 
cent  works  which  report  local  current  density  data  using  seg¬ 
mented  current  collectors  [1],  typically  the  experimental  data 
are  available  only  for  global  output  quantities  such  as  polar¬ 
ization  curves,  overall  water  balances,  etc.  Furthermore,  these 
data  are  actually  a  set  of  measurements  over  a  range  of  oper¬ 
ating  conditions  as  a  parameter,  such  as  the  average  current 
density,  cell  potential,  relative  humidity,  etc.  is  varied.  There¬ 
fore,  the  comparison  between  model  and  experiment  can  only 
be  made  globally  using  data  that  represent  an  average  perfor¬ 
mance  of  the  fuel  cell.  This  point  is  particularly  significant 
in  the  case  of  fuel  cell  stacks,  where  the  local  current,  mass 
fluxes,  temperature,  etc.  may  vary  significantly  from  the  av¬ 
erage  measured  data. 

For  the  case  of  polymer  electrolyte  membrane  fuel  cells 
(PEMFCs),  the  mathematical  model  is  chosen  to  represent 
various  reaction  and  transport  processes  within  a  cell,  such  as 
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convection-diffusion  transport,  electrochemical  reaction  and 
membrane  water  and  proton  transport.  The  solution  variables 
can  include  fluid  velocity  and  pressure,  species  composition, 
temperature,  scalar  potentials,  etc.  while  the  parameters  may 
include  the  membrane  conductivity,  reference  exchange  cur¬ 
rent  densities,  electro-osmotic  drag  coefficient,  etc.  An  ad¬ 
ditional  motivation  for  parameter  estimation  is  the  lack  of 
understanding  of  fundamental  transport  processes  such  as  mi¬ 
gration  of  protons  and  liquid  water  through  the  membrane, 
which  can  lead  to  semi-empirical  models  with  unknown, 
sometimes  nonphysical  parameters,  which  must  somehow  be 
estimated. 

Many  mathematical  models  for  PEMFCs  have  been  de¬ 
scribed  in  the  literature,  ranging  from  one-dimensional  mod¬ 
els  in  the  direction  of  current  or  gas  flow  [1-5]  to  two-  and 
three-dimensional  models  [6-10].  However,  most  of  these 
authors  do  not  address  the  question  of  parameter  estimation 
using  a  systematic  approach.  Thampan  et  al.  [11]  and  Fim- 
rite  [12]  performed  parameter  estimation  for  membrane  con¬ 
ductivity  submodels  using  curve  fitting  techniques.  Suares 
and  Hoo  [13]  estimated  PEMFC  model  parameters  for  the 
model  of  Nguyen  and  White  [4],  such  as  exchange  cur¬ 
rent  density  using  an  optimization  approach.  Berg  et  al. 
[1]  estimate  several  parameters  using  a  one-dimensional 
PEMFC  model  based  on  [4].  Recently,  Guo  et  al.  [14] 
have  fitted  cathode  catalyst  layer  parameters  such  as  porosi¬ 
ties,  reference  current  densities  and  effective  diffusion  co¬ 
efficients  using  a  one-dimensional  cathode  catalyst  layer 
model. 

In  this  work,  we  demonstrate  the  use  of  parameter  estima¬ 
tion  for  a  one-dimensional  model  for  a  PEMFC  that  accounts 
for  conduction  of  protons  and  electrons,  diffusion  of  oxygen 
and  the  electrochemical  reactions.  We  present  results  on  esti¬ 
mation  of  the  effective  membrane  conductivity,  effective  vol¬ 
umetric  exchange  current  densities  in  the  cathode  and  anode 
catalyst  layers,  and  effective  oxygen  diffusion  coefficients  in 
the  gas  diffusion  and  catalyst  layers,  paying  special  attention 
to  the  sensitivity  coefficients,  which  indicate  the  significance 
of  a  parameter  at  each  data  point  on  the  experimental  po¬ 
larization  curve.  The  algorithm  for  parameter  estimation  is 
based  on  least  squares  minimization  with  constraints  arising 
from  the  solution  of  a  system  of  nonlinear  PDEs,  and  is  suf¬ 
ficiently  general  that  it  can  be  extended  to  multidimensional, 
transient  models,  or  to  include  additional  physical  phenom¬ 


ena  in  the  mathematical  model.  In  addition,  we  study  the 
effects  of  mesh  refinement  on  the  accuracy  of  the  parameter 
estimation,  including  the  use  of  locally  refined  meshes  for 
greater  efficiency  and  reliability. 

The  remainder  of  this  work  is  organized  as  follows.  We 
first  describe  our  one-dimensional  transport  model  for  a 
PEMFC.  Then  we  provide  an  algorithm  for  parameter  estima¬ 
tion  with  PDE  constraints  based  on  a  nonlinear  least  squares 
approach.  Next,  we  present  numerical  results  for  estimating 
parameters  using  conventional  PEMFC  dimensions  and  data 
for  regions  with  no  mass  transport  limitations,  as  well  as  for 
an  air-breathing  PEMFC  that  is  mass  transport  limited.  We 
conclude  with  an  example  illustrating  the  effectiveness  of 
using  adaptive  mesh  refinement. 


2.  One-dimensional  PEMFC  model  equations 

In  this  section,  we  define  a  simple  one-dimensional  model 
for  reaction  and  transport  within  a  PEMFC.  Since  the  prin¬ 
cipal  axis  of  transport  of  current  is  in  the  direction  normal 
to  the  membrane,  we  choose  the  v-axis  in  this  direction,  ori¬ 
ented  from  anode  to  cathode.  A  schematic  of  the  PEMFC 
layers  and  dependent  variables  is  shown  in  Fig.  1. 

Conduction  of  electrons  through  the  catalyst  layers  and  gas 
diffusion  layers  (GDFs)  is  modeled  using  a  solid  electrical 
potential  0e  and  conduction  of  the  protons  in  the  membrane 
and  catalyst  layers  is  modeled  using  an  ionic  potential  0m. 
In  the  cathode  side  GDF  and  catalyst  layer,  we  model  the 
transport  of  oxygen  using  the  oxygen  mole  fraction  Xo2. 
Conservation  of  charge  and  mass,  combined  with  Ohm’s  law 
and  Fick’s  law,  yields  the  nonlinear  boundary  value  problem: 
find  {(pe,  cf>m,  Xo2}  such  that 
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Fig.  1.  Schematic  of  one-dimensional  PEMFC  model. 
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r]  =  <fiQ  —  0m,  overpotential  in  catalyst  layers.  (Id) 

Here  crQ  is  the  electrical  conductivity,  crm  the  ionic  conductiv¬ 
ity,  c  the  molar  concentration  and  D  is  the  molecular  diffu- 
sivity.  The  nonlinear  reaction  kinetics  are  given  by  a  Butler- 
Volmer  expression 

j  =  is  {exp  (pnfil)  -  exp  ((0  -  1  )nfrj)} ,  (2) 

where  js  is  the  exchange  current  density  with  s  =  a,  c  for 
anode  and  cathode,  respectively,  ft  the  transfer  coefficient,  n 
the  number  of  electrons  transferred  and  /  =  F/RT  is  a  factor 
depending  on  the  Faraday  constant  F,  ideal  gas  constant  R 
and  temperature  T  [15].  A  common  assumption  [10]  is  that 
P  =  1/2,  in  which  case  we  can  define  a  =  /3n  =  —  (ft  —  \)n 
and  then 


j  =  is  {exp  (afrj)  -  exp  (-afy)}  =  js 2  sin  h(afy).  (3) 

In  general,  the  transport  of  hydrogen  in  the  anode  GDL  and 
consumption  in  the  anode  catalyst  layer  should  be  included  in 
the  model.  However,  in  the  context  of  fitting  parameters  using 
polarization  curve  data,  the  hydrogen  mass  transport  limita¬ 
tions  and  resulting  overpotential  losses  are  likely  to  be  dom¬ 
inated  by  ohmic,  cathodic  overpotential  and  oxygen  limiting 
losses  [6].  It  is  therefore  reasonable  to  assume  that  hydrogen 
concentration  is  nearly  constant,  while  oxygen  concentration 
is  expected  to  vary  more  widely.  We  make  the  assumption 
that  the  exchange  current  density  is  constant  in  the  anode  and 
first-order  in  the  cathode 
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The  boundary  conditions  are  specified  as  follows: 


anode  GDL  boundary  with  the  anode  gas  channel:  0e  =  0; 
GDL/catalyst  layer:  zero  ionic  current  or  — crm^0m  =  0; 
catalyst  layer/membrane:  zero  electrical  current  or 

— <Te^</>e  =  0; 

membrane/cathode  catalyst  layer:  zero  oxygen  flux  or 

-cDfx  02  =  0; 

cathode  GDL  boundary  with  cathode  gas  channel:  total  po¬ 
tential  drop  0e  =  —  dL  and  reference  oxygen  mole  fraction 

V  _  vambient 

A02  —  aq2 


We  now  define  a  subset  of  parameters  on  which  to  apply 
the  parameter  estimation  procedure  of  Section  3.  These  are 
the  membrane  ionic  conductivity  amem,  the  cathode  reference 
exchange  current  density  j‘Jef,  the  anode  reference  exchange 
current  density  j^ef,  and  the  oxygen  diffusion  coefficients 
Dgdi  and  Dcat  in  the  cathode  side  GDL  and  catalyst  layer, 
respectively.  We  will  assume  that  all  other  parameters  and 
dimensions  are  known,  and  will  seek  to  estimate  only  these 
five  parameters. 


3.  Algorithm  for  nonlinear  least  squares  with  PDE 
constraints 

In  this  section,  we  briefly  outline  the  nonlinear  least 
squares  (NLS)  approach,  introduce  the  concept  of  implicit 
least  squares  through  PDE  constraints,  and  give  an  algorithm 
for  NLS  with  PDE  constraints. 

The  Lagrangian  approach  to  constrained  optimization  is 
well  known  (cf.  [16,17]).  As  part  of  this  work,  an  algorithm 
was  developed  by  the  authors  based  on  extending  the  La¬ 
grangian  approach.  This  allows  for  elimination  of  the  La¬ 
grange  multipliers,  avoiding  solution  of  a  large  coupled  sys¬ 
tem  for  the  parameters,  solutions,  and  Lagrange  multipliers. 
In  exchange,  we  must  solve  additional  linear  equations  for 
the  sensitivity  functions,  but  can  sequentially  solve  a  number 
of  decoupled  systems  in  order  to  compute  each  parameter 
update. 

Suppose  that  we  have  M  data  pairs  {(/;,  V/)}^1  and  a  fit¬ 
ting  function  I  =  /(A;  V)  with  N  parameters  X  =  {Ay}^. 
For  example,  these  may  be  taken  as  the  current  density  I  and 
cell  voltage  V  from  polarization  curve  data,  as  will  be  done  in 
Section  4.  We  would  like  to  choose  the  parameters  A  in  order 
to  approximate  the  data  at  each  data  pair  as: 

I(KVi)^Ii.  (5) 

A  common  way  to  do  this  is  to  minimize  a  least  squares 
functional 

1  1  M 

M{X)  =  -\\m\\2  =  -  ^  \mi(X)\2, 

i= 1 

mm  =  Vi)  -  Ii),  (6) 


For  conventional  PEMFCs  with  flow  fields,  x^bient  will  vary 
along  the  flow  channels,  and  can  be  interpreted  as  an  average 
value,  while  for  ambient  air-breathing  PEMFCs,  we  can  de¬ 
fine  xg£bient  using  the  ambient  oxygen  composition.  In  the 
fitting  procedure  later,  we  only  consider  the  oxygen  equation 
for  ambient  air-breathing  PEMFCs. 

We  specify  a  given  potential  drop  dV,  solve  the  model 
equations,  and  then  calculate  the  current  density  by  com¬ 
puting  the  electrical  current  /  =  —  cre^(/>e  evaluated  at  the 
cathode  GDL  boundary.  The  cell  voltage  is  then  computed  to 
be  V  =  Voc  —  dV,  where  Voc  is  the  experimentally  measured 
open  circuit  potential. 


using  a  NLS  algorithm,  where  the  mf  s  are  called  the  misfits. 
The  constants  coj  can  be  chosen  to  produce  a  weighted  data 
fit.  For  simplicity  we  take  coi  =  1  throughout. 

At  a  local  minimum  we  must  have  VA4  =  (Vm)T  m  =  0, 
or 


which  defines  the  set  of  solutions  of  the  NLS  problem.  In 
order  to  solve  (7),  we  can  apply  Newton-Raphson  iteration. 
The  system  that  we  solve  for  the  update  SX  at  each  iteration 
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takes  the  form  J  9  A  =  r,  where 


subject  to  the  constraints 


dm  i 

9V 


M 

E 

i=l 


dmi  dm , 


dXj  dXk 


+  nti 


d2mi 

dXjdXk 


If  we  assume  each  mz-  is  small,  we  can  drop  the  second  term 
and  solve  using  the  Gauss-Newton  method,  which  takes  the 
form  M  9  A  =  r,  where 


E 


dmi  9m; 
dXj  dXk 


One  advantage  to  this  approach  is  that  we  only  need  to  calcu- 
late  the  first-order  partial  derivatives  and  no  second-order 

derivatives. 

The  partial  derivatives  are  often  called  the  sensitiv¬ 
ity  coefficients,  because  they  measure  the  sensitivity  of  the 
ith  misfit  to  perturbations  in  the  jth  parameter.  They  are  also 
important  to  the  question  of  solvability  of  the  NLS  problem, 
which  is  also  called  the  problem  of  identifiability  (see  for  in¬ 
stance  [16]).  Invertibility  of  the  matrix  M  is  equivalent  to  the 
vectors  a-7  =  ,  . . . ,  G  Mm  being  linearly  indepen¬ 

dent.  If  the  vectors  are  nearly  linearly  dependent,  then  the 
NLS  problem  can  be  difficult  to  solve  and  may  even  fail  to 
converge.  Thus,  before  proceeding  with  any  NLS  problem, 
it  is  prudent  to  identify  any  linear  dependence  in  the  vectors 
aJ . 

The  NLS  approach  described  above  is  adequate  when  we 
have  an  explicit  expression  of  the  form  I  =  /(A;  V).  However, 
often  we  are  interested  in  solving  stationary  or  transient  PDEs 
with  specified  input  value  V  and  postprocessing  the  solution 
to  obtain  the  output  value  /,  or  vice  versa.  The  parameters 
A  are  then  the  parameters  of  the  PDE,  and  the  “function” 
I  =  /(A;  V )  is  computed  by  solving  the  PDE. 

In  this  case,  we  must  modify  the  NLS  approach  to  in¬ 
clude  the  PDE  constraints  and  to  handle  the  implicitly  de¬ 
fined  output  function.  If  we  denote  the  differential  operator 
that  models  the  the  PEMFC  by  R  =  R(u,  A;  V),  then  solu¬ 
tions  u  =  u( A;  V)  are  defined  by  solving: 


R(u,  A;V)  =  0, 


(10) 


for  given  values  of  A  and  V.  The  output  /,  is  calculated  as  a 
functional  of  the  solution,  in  the  form  I  =  I(u).  We  denote 
the  residuals  for  each  value  VJ  by  Rl(u,  A)  =  R(u,  A;  Vf)  with 
corresponding  solutions  mz  =  mz( A). 

The  constrained  parameter  estimation  problem  can  then 
be  expressed  as  follows:  find  u  and  A  that  minimizes  the 
functional 


Ai(u)  = 


1 

2 


M 

E  N<(wi)l2, 

i=  1 


mt(ui)  =  I(m)  -  It, 


(11) 


Ri(ui,  A)  =  0. 


(12) 


The  main  difference  is  that  the  misfit  functional  is  now  a 
function  of  the  solution  variable  u  rather  than  the  parameter 
A.  The  dependence  of  u  on  A  is  implicit  in  the  constraints. 
We  will  assume  for  simplicity  that /is  linear,  so  that  I\u)v  = 

m. 

A  standard  approach  to  solving  the  constrained  parameter 
estimation  problem  is  to  introduce  Lagrange  multipliers,  one 
for  each  solution  mz,  and  solve  a  large  coupled  system  [17].  It 
turns  out  that  we  can  eliminate  the  Lagrange  multipliers  and 
solve  the  constrained  parameter  estimation  problem  using  a 
slightly  modified  form  of  the  Gauss-Newton  iteration  in  (9), 
provided  that  we  can  calculate  the  partial  derivatives  and 

the  components  m;.  Using  the  chain  rule,  we  can  calculate 


dmi 

dXj 


dmi  dui 
dui  9  A  j 


(13) 


Thus,  if  we  could  calculate  the  partial  derivatives  Jp,  we 

could  evaluate  .  To  see  how  to  obtain  these  functions,  we 
differentiate  each  PDE  w.r.t.  Xj  to  obtain: 


3R'  dui  3Rl 

- -  + - =0. 

du  i  3 Xj  3 Xj 

Thus,  by  solving  the  auxiliary  equations: 


(14) 


3  Rl  3  Rl 

XX)  i  j  —  -  , 

dui  9  A  j 


(15) 


we  obtain  Wij  =  and  thus  =  I(iVij). 

We  now  define  an  algorithm  for  parameter  estimation  with 
PDE  constraints  shown  in  the  box  below. 

Because  the  system  in  step  (2)  involves  summation  over 
the  data  index  /,  we  can  perform  each  iteration  of  the  algo¬ 
rithm,  which  consists  of  steps  (1-3),  by  looping  over  i  and 
then  j.  Practically  speaking,  this  means  that  we  need  only 
form  and  store  the  Jacobian  matrix  associated  with  the  dif- 

r) 

ferential  operator  ^  once  for  each  index  i  in  step  (1)  of  each 
iteration.  The  linear  systems  can  be  solved  by  an  LU  factor¬ 
ization  of  the  Jacobian  matrix  or  else  an  iterative  solver  can  be 
applied.  The  finite  dimensional  approximate  solution  vectors 
associated  with  8ui  and  are  stored  for  the  entire  iteration 
for  use  in  step  (3).  Thus,  the  cost  of  each  iteration  of  the  algo¬ 
rithm  is  the  same  as  solving  N(M  +1)  Jacobian  systems.  The 
additional  solution  of  the  Gauss-Newton  system  in  step  (2) 
is  then  negligible.  The  result  is  the  simultaneous  calculation 
of  the  optimal  parameters  A  j,  the  corresponding  solutions  Ui, 

r\ 

the  misfit  values  mz,  and  the  sensitivity  coefficients  . 

This  algorithm  can  be  implemented  in  any  computer  code 
which  allows  access  to  the  Jacobian  matrix  and  residual  vec- 
tor  and  the  linear  algebra  solver.  If  the  partial  derivatives  ^ 

cannot  be  calculated  explicitly,  we  may  approximate  them 


B.  Carnes,  N.  Djilali  /  Journal  of  Power  Sources  144  (2005)  83-93 


87 


Algorithm  for  NLS  with  PDE  constraints 

Initialize  mz-  for  initial  A  by  solving  Rl(ui,  A)  =  0. 
Compute  updates  8ut  and  8k  j  as  follows: 
for  n  =  1,  2,  ...  do 


(1)  Compute  mi  =  I(uf)  —  /;.  Solve  for  initial  8ui  and 


functions  w 


ij 


dRl 

dui 

dR1 

du; 


8ui  =  —  Rl  (Newton-Raphson  iteration).  (16a) 


Wij  =  - 


dR1 
dk ; 


(16b) 


(2)  Solve  M  8k  =  r  using  modified  Gauss-Newton  iter¬ 
ation 


M 


fj  =  -  ffinii  +  I(8ui)}I(wij), 


i= 1 


M 


Mjk  —  ^ 


i=l 


(3)  Update  8ui  using  correction  for  A 


N 


8u[+  —  ^  ^  WijSXj. 
7=1 


(17a) 


(17b) 


(18) 


Update  Ui+  =  8ut  and  kj+  =  8k  j,  possibly  with  step 
size  control. 


end  for 


using  the  finite  difference  approximation: 

«  -[R\ui,  A  +  eej)  -  R\ui,  A)],  (19) 

dkj  £ 

where  ej  is  the  jth  coordinate  vector  and  s  1 .  If  the  code 
does  not  permit  access  to  the  Jacobian  matrix,  then  we  may 
directly  approximate  the  partial  derivatives  using: 

3wi  ■  1 

—f  -  \nii(ui(X  +  eej))  -  mj(ui(X))] ,  (20) 

where  w;(A)  is  the  solution  calculated  at  the  parameter  value 
A. 

We  now  apply  the  parameter  estimation  algorithm  to  our 
PEMFC  model  equations. 

4.  Application  to  the  PEMFC  model  equations 

The  solution  to  the  differential  equations  is  approximated 
at  each  step  of  the  parameter  estimation  algorithm  by  a  stan¬ 
dard  Galerkin  finite  element  method  using  piecewise  lin¬ 
ear  basis  functions.  We  begin  with  an  initial  coarse  grid 


and  apply  the  iterative  algorithm,  which  usually  converges 
in  about  four  to  eight  iterations,  to  solve  for  approximate 
solutions  Uih  and  parameters  A^.  Then  we  refine  the  grid 
and  project  the  current  approximate  solutions  onto  the  fine 
grid.  Using  the  projected  solutions  as  the  initial  guess  we 
apply  the  iterative  algorithm  to  obtain  new  approximate  so¬ 
lutions  and  parameters.  After  sufficient  levels  of  mesh  refine¬ 
ment,  the  parameters  tend  to  converge  to  a  mesh-independent 
value,  which  we  take  as  the  final  estimated  parameter 
value. 

The  mesh  refinement  can  be  uniform,  meaning  that  all  cells 
are  refined,  or  adaptive,  meaning  that  we  selectively  refine 
cells  based  on  a  local  error  indicator  function.  Our  computer 
code  reports  the  values  of  the  parameters  calculated  at  each 
iteration  for  each  mesh,  which  reveals  the  convergence  of  the 
parameter  estimation  algorithm  on  each  individual  mesh,  as 
well  as  the  convergence  of  the  finite  element  approximation 
to  the  parameters. 

We  consider  data  from  two  different  PEMFCs,  which 
are  shown  in  Fig.  2.  The  first  data  set  is  for  a  conven¬ 
tional  PEMFC  design  using  polarization  data  reported  by 
Ticianelli  et  al.  [18],  while  the  second  is  an  air-breathing 
PEMFC  yielding  lower  current  densities.  We  also  plot  our 
best  fitted  polarization  curves,  corresponding  to  the  simulta¬ 
neous  fitting  of  (Tmem  and  jc  for  the  conventional  PEMFC 
and  to  (Tmem,  T>cat,  and  Dg di  for  the  low  current  density 
PEMFC. 

In  the  parameter  estimation  procedure,  we  can  choose  to 
fit  parameters  individually,  or  simultaneously.  While  the  for¬ 
mer  approach  is  numerically  more  reliable,  it  may  not  rep¬ 
resent  the  coupling  between  the  parameters.  Therefore,  the 
latter  approach  is  preferred  based  on  physical  grounds.  We 
choose  to  perform  both  approaches,  in  order  to  demonstrate 
the  variation  in  the  fit  of  the  parameters.  When  a  great  deal  of 
variation  of  a  fitted  parameter  is  seen  using  different  combi¬ 
nations,  this  can  be  an  indication  that  there  is  strong  coupling 
between  the  parameters.  In  this  case  the  data  can  be  fit  by 
multiple  sets  of  parameters,  and  the  model  may  be  need  to  be 
revised. 

4.1.  Conventional  PEMFC 

We  first  fit  our  one-dimensional  PEMFC  model  to  the  con¬ 
ventional  PEMFC  polarization  curve. 

4.1.1.  Initial  parameters 

We  take  the  physical  dimensions  of  each  layer  and  ini¬ 
tial  transport  and  reaction  properties  from  Berning  et  al.  [6], 
which  we  summarize  in  Table  1.  The  experimental  data  are 
from  Ticianelli  et  al.  [18]  and  takes  the  form  of  a  curve  fit 

v=  Vo-Mog(2j  _RIt  (21) 

with  parameter  values  Vo  =  0.935  V,  b  =  0.065  V,  7o  = 
1mA  cm-2,  and  R  =  0.39  cm2  S-1.  The  electrode  is  pro¬ 
duced  with  a  catalyst  loading  of  0.35  mg  Pt  cm-2  with  4% 


88 


B.  Carnes,  N.  Djilali  /  Journal  of  Power  Sources  144  (2005)  83-93 


Current  density  [mA/cmA2] 

Fig.  2.  Experimental  and  fitted  polarization  curves  for  conventional  and  low 
current  PEMFC  designs. 


Nafion,  the  membrane  is  Nafion  117,  and  Voc  =  0.935.  The 
PEMFC  is  operating  on  humidified  hydrogen/air  at  pressures 
of  3/5  atm  at  80  °C.  Because  the  data  are  fitted  to  lower  current 
density  data,  mass  transport  limitations  are  not  significant. 
Therefore,  we  do  not  include  the  oxygen  transport  equation, 
assuming  that  the  oxygen  concentration  is  at  the  reference 
value.  We  seek  to  fit  the  parameters  <rme m,  and  ja  (drop¬ 
ping  the  superscript  on  jJef). 

For  fully  humidified  Nafion  117,  we  expect  <7mem  ~ 
8  S  m_  1  from  experimental  measurements  of  the  protonic 
conductivity  of  Nafion.  We  choose  the  initial  value  of 
(Tmem  =  6Sm_1  following  [6].  The  exchange  current  den¬ 
sities  are  typically  available  as  surface  reaction  rates;  for 
example,  in  [6]  they  use  values  of  ja  =  6.0e+3  Am-2  and 
jc  =  4.4e— 3  Am-2.  In  order  to  derive  effective  volumetric 
exchange  current  densities,  we  must  multiply  by  a  correc¬ 
tion  factor  that  has  units  of  1  m-1.  Ideally  this  factor  should 
account  for  the  surface:  volume  ratio,  the  porosity  of  the  cat¬ 
alyst  layer,  as  well  as  the  amount  of  Nafion  and  Pt  applied 
when  the  catalyst  layer  was  formed.  We  take  our  initial  guess 
for  the  effective  volumetric  exchange  current  densities  as 
ja  =  le+9  Am-3  and  jc  =  2e+6Am-3. 

Table  1 


Initial  physical  transport  and  reaction  parameters  for  conventional  PEMFC 


Anode 

Membrane 

Cathode 

Gas  diff. 
layer 

Cat. 

layer 

Cat. 

layer 

Gas  diff. 
layer 

Thickness,  l  (fxm) 

260 

28 

230 

28 

260 

o'e  (S  m-1) 

6000 

6000 

0 

6000 

6000 

(jm  (S  m-1) 

0 

0.6 

6 

0.6 

0 

Js  (Am-3) 

le+9 

2e+6 

a 

0.5 

1.0 

We  present  in  Fig.  3  the  dimensionless  sensitivity  coef- 
ficients  for  all  three  parameters,  calculated  at  the  initial 

values.  Because  of  the  near  linear  dependence  of  the  sen¬ 
sitivity  coefficients  for  txmem  and  ja,  we  expect  that  these 
parameters  may  be  coupled. 

Another  interesting  observation  to  make  about  Fig.  3 
is  the  variation  in  sensitivity  as  the  average  current  den¬ 
sity  I  is  varied.  While  the  sensitivity  coefficients  for  crm 
and  ja  remain  fairly  constant,  the  sensitivity  coefficients 
for  jc  are  largest  for  small  /,  and  decrease  exponentially 
as  I  increases.  This  can  be  explained  by  the  difference  in 
the  type  of  voltage  losses.  For  the  membrane,  the  ohmic 
loss  is  a  linear  function  of  the  current  with  constant  sen¬ 
sitivity  to  the  conductivity,  while  for  the  cathode  electro¬ 
chemical  reaction,  the  loss  is  a  logarithmic  function  of  the 
current  with  decreasing  sensitivity  to  the  reaction  rate  as 
I  increases.  By  linearizing  the  anode  reaction  rate,  it  can 
be  shown  that  the  anode  loss  is  also  a  linear  function  of 
the  current,  and  thus,  behaves  as  the  membrane  conductiv¬ 
ity. 

In  many  studies,  the  catalyst  layers  are  modeled  as  thin 
interfaces  or  using  a  constant  overpotential.  In  Fig.  4  we 
plot  the  overpotential  rj  in  the  anode  and  cathode  cata¬ 
lyst  layers  for  the  initial  parameter  values  and  varying  cur¬ 
rent  densities  in  order  to  illustrate  the  effect  of  solving  a 
one-dimensional  model  for  behavior  that  cannot  be  mod¬ 
eled  using  a  reduced  model.  The  maximal  absolute  val¬ 
ues  of  rj  occurs  along  the  membrane-catalyst  layer  inter¬ 
faces,  where  the  rate  function  has  its  largest  absolute  value. 
The  larger  rate  constant  at  the  anode  results  in  faster  ki¬ 
netics  and  a  smaller  reaction  region  near  the  interface.  We 
do  not  plot  the  solutions  elsewhere,  since  0m  varies  lin¬ 
early  within  the  membrane,  and  0e  is  nearly  constant  in  the 
GDFs. 

4.1.2.  Estimated  parameters 

We  now  estimate  the  parameters  using  our  parameter  es¬ 
timation  algorithm.  We  first  compare  the  results  obtained  by 


Fig.  3.  Initial  dimensionless  sensitivity  coefficients  for  conventional 
PEMFC. 
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Fig.  4.  Initial  overpotential  in  anode  (left)  and  cathode  (right).  I  is  given  in  units  of  mA  cm 


Table  2 


Estimated  parameters  on  fine  grids  for  conventional  PEMFC  using  different 
combinations  of  fitted  parameters 


^mem  (S  m  ) 

jc  (Am  3) 

ja  (A m  3) 

Initial 

6 

2e+6 

le+9 

Fitted 

7.67 

— 

— 

— 

9.79e+5 

— 

(7.67) 

— 

1.43e+9 

8.29 

2.50e+6 

— 

9.24 

— 

9.60e+8 

(7.67) 

2.66e+6 

9.02e+8 

Values  in  parenthesis  were  fixed  but  chosen  different  from  the  initial  guesses. 


coefficient  is  small,  the  term  in  (7)  multiplied  by  the  coeffi¬ 
cient  will  also  be  small,  and  thus  can  be  neglected.  Therefore 
we  choose  to  estimate  jc  using  current  density  only  in  the 
range  of  I  <  300  mA  cm-2,  where  the  sensitivity  is  largest. 
Inclusion  of  jc  is  thus  key  to  fitting  low  values  of  /,  whereas 
inclusion  of  <7mem  is  required  for  the  entire  range  of  I.  This 
is  consistent  with  the  dominance  of  activation  losses  at  low 
current  densities  and  of  ohmic  losses  at  higher  current  den¬ 
sities,  and  explains  why  the  best  fittings  are  obtained  using 
the  combinations  {<Tmem,  jc)  and  { jc,  ja}  (with  the  fitted  <7mem 
value). 


estimating  each  parameter  separately,  and  then  using  combi¬ 
nations  of  parameters.  The  converged  values  are  reported  in 
Table  2. 

Convergence  of  the  estimated  values  for  amem  using  uni¬ 
form  mesh  refinement  is  shown  in  Fig.  5.  We  can  see  con¬ 
vergence  of  the  parameter  estimation  algorithm  on  each  grid 
to  a  value  Ajh,  as  well  as  convergence  of  the  values  Ajh  as 
the  mesh  is  refined.  When  we  fit  <7mem  or  both  <7mem,  jc,  the 
value  is  close  to  the  fully  humidified  conductivity  of  8  S  m-1 
for  Nafion  117.  However,  fitting  with  both  amem,  ja  yields  a 
higher  value  of  crmem  hy  fitting  a  lower  value  of  ja. 

We  saw  in  Fig.  3  that  the  sensitivity  of  the  current  density  I 
to  jc  decreases  dramatically  as  /  increases.  When  a  sensitivity 


Table  3 


Initial  physical  transport  and  reaction  parameters  for  a  low  current  PEMFC 


Anode 

Membrane 

Cathode 

Gas  diff. 
layer 

Cat. 

layer 

Cat. 

layer 

Gas  diff. 
layer 

Thickness,  /  (|xm) 

130 

20 

60 

20 

130 

cre  (S  m-1) 

100 

100 

0 

100 

100 

crm  (S  m-1) 

0 

0.2 

2 

0.2 

0 

D0l  (m2  s—  1 ) 

0 

0 

0 

3e— 7 

3e— 6 

js  (Am-3) 

le+9 

2e+6 

a 

0.5 

1.5 

4.2.  Low  current  density  PEMFC 

We  now  fit  our  one-dimensional  PEMFC  model  to  the  low 
current  density  PEMFC  polarization  curve. 


Fig.  5.  Convergence  of  amein  for  various  parameter  combinations. 


90 


B.  Carnes,  N.  Djilali  /  Journal  of  Power  Sources  144  (2005)  83-93 


Fig.  6.  Initial  dimensionless  sensitivity  coefficients  for  a  low  current 
PEMFC. 

4.2.1.  Initial  parameters 

We  adjust  the  dimensions  and  parameters  for  an  air- 
breathing  PEMFC  which  exhibits  mass  transport  limitations 
at  relatively  low  current  densities  as  shown  in  Fig.  2.  This 
cell  uses  a  thinner  membrane  (Nation  112)  and  thinner  lay¬ 
ers,  which  we  summarize  in  Table  3.  The  temperature  is  now 
specified  as  60  °C,  V0c  =  0.912,  and  the  ambient  oxygen 
mole  fraction  is  specified  as  0.21.  Because  of  the  reduced 
slope  of  the  polarization  curve,  we  begin  with  a  lower  value 
for  the  membrane  conductivity.  The  oxygen  diffusion  coef¬ 
ficient  in  the  GDF  is  reduced  by  an  order  of  magnitude  to 
account  for  the  porosity  and  tortuosity  of  the  diffusion  layer; 
in  the  catalyst  layer  this  parameter  is  reduced  by  another  or¬ 
der  of  magnitude  to  model  additional  resistance  at  reaction 
sites  (see  Fig.  6  for  initial  sensitivites). 

4.2.2.  Estimated  parameters 

We  now  apply  the  parameter  estimation  to  the  data  for 
the  low  current  PEMFC  design.  The  converged  estimated 
parameters  are  shown  in  Table  4. 


The  fitting  of  the  parameters  <7mem,  Ja>  and  jc  is  per¬ 
formed  using  only  data  with  I  <  200  mA  cm-2  and  is  sim¬ 
ilar  to  the  case  of  the  conventional  fuel  cell  (see  first 
section  of  fitted  parameters  in  Table  4).  However,  the 
fitted  membrane  conductivity  is  now  closer  to  ISm-1, 
which  corresponds  to  partially  humidified  operating  con¬ 
ditions.  In  addition,  the  fitted  cathode  rate  jc  is  an  or¬ 
der  of  magnitude  lower  than  the  initial  value,  which  is 
caused  by  using  the  effective  value  of  2e—  ISm-1  for 
the  effective  membrane  conductivity  in  the  catalyst  lay¬ 
ers. 

We  next  fit  the  effective  oxygen  diffusion  coefficients  Dg&\ 
and  Dcat  (second  section  of  fitted  parameters  in  Table  4)  for 
the  large  current  data  (/  >  200  mA  cm-2)  using  the  fitted 
ffmem?  ja,  and  jc  values.  The  data  can  be  fitted  by  reducing 
either  Dg di  or  Dcat  by  an  order  of  magnitude,  representing 
an  additional  mass  transport  resistance,  such  as  liquid  wa¬ 
ter  flooding.  We  can  also  fit  the  diffusion  coefficients,  indi¬ 
vidually  or  simultaneously,  along  with  crmem,  over  the  en¬ 
tire  polarization  curve,  using  the  fitted  values  for  ja  and  jc 
(third  section  of  fitted  parameters  in  Table  4).  We  plot  in  Fig. 
7  the  oxygen  profiles  in  the  cathode  catalyst  layer  for  dif¬ 
ferent  current  densities  and  parameter  fittings.  Fitting  Dgc \\ 
results  in  a  nearly  constant,  low  oxygen  concentration  in 
the  catalyst  layer,  while  fitting  Dcat  results  in  a  steep  gra¬ 
dient  in  oxygen  concentration  in  the  catalyst  layer.  Simulta¬ 
neous  fitting  lies  somewhere  in  between,  and  distributes  the 
mass  transfer  resistance  between  the  GDF  and  the  catalyst 
layer. 

While  simultaneous  fitting  provides  the  best  possible  data 
fit  over  the  entire  curve,  all  three  fits  are  similarly  good.  The 
problem  lies  not  in  the  parameter  estimation  procedure,  but 
in  the  choice  of  the  mathematical  model.  In  order  to  fit  the 
drop  in  current  density  of  the  data  using  the  oxygen  equa¬ 
tion,  we  must  add  a  certain  amount  of  resistance  to  mass 
transport,  in  the  form  of  reduced  diffusion  coefficients  Dg&\ 
and  Dcat.  The  fairly  small  diffusion  coefficients  required  to 
obtain  a  good  fit  (two  to  three  orders  of  magnitude  smaller 
than  for  pure  gas  mixtures),  suggest  that  a  macrohomoge- 
neous  mass  transport  resistance  model  may  be  inadequate  to 


Table  4 


Estimated  parameters  for  fine  grids  for  low  current  PEMFC  using  different  combinations  of  fitted  parameters 


^mem  (S  m  ) 

jc  (Am  3) 

ja  (A m  3) 

Dgdi(m2s  l) 

Dcat(m2s  l) 

Initial 

2 

2e+6 

le+9 

3e— 6 

3e— 7 

Fitted 

0.715 

— 

— 

— 

— 

— 

3.71e+5 

— 

— 

— 

1.02 

3.41e+5 

— 

— 

— 

— 

— 

8.27e+7 

— 

— 

(1.02) 

(3.41e+5) 

1.03e+9 

— 

— 

Fitted 

(1-02) 

(3.41e+5) 

(1.03e+9) 

1.66e— 7 

— 

(1.02) 

(3.41e+5) 

(1.03e+9) 

— 

1.89e— 8 

Fitted 

1.10 

(3.41e+5) 

(1.03e+9) 

1.63e— 7 

— 

1.16 

(3.41e+5) 

(1.03e+9) 

— 

1.76e— 8 

1.15 

(3.41e+5) 

(1.03e+9) 

5.27e— 7 

2.38e— 8 

Values  in  parenthesis  were  fixed  but  chosen  different  from  the  initial  guesses. 


B.  Carnes,  N.  Djilali  /  Journal  of  Power  Sources  144  (2005)  83-93 


91 


Fig.  7.  Oxygen  mole  fraction  profiles  in  the  cathode  catalyst  layers  obtained  by  fitting  (Jmem,  Dc at  and  Dg(ji  for  I  =  3 17  mA  cm  2  (left)  and  I  =  356  mA  cm 
(right). 


describe  the  current  losses,  and  an  improved  model  should 
be  sought. 

Finally,  we  note  that  for  large  current  densities  where  the 
model  predicts  severe  oxygen  depletion,  the  reaction  rate  can 
even  develop  an  interior  maximum  within  the  catalyst  layer, 
instead  of  along  the  membrane/catalyst  layer  interface,  as 
seen  in  Fig.  8. 

4.3.  Adaptive  mesh  refinement 

In  this  section,  we  demonstrate  the  use  of  adaptive  mesh 
refinement  (AMR)  to  increase  the  computational  efficiency  of 
the  parameter  estimation  algorithm.  Here  instead  of  refining 
all  elements  of  the  mesh  to  obtain  the  next  finer  mesh,  we 
choose  only  a  subset  of  elements  for  refinement.  To  do  this 
we  need  an  error  indicator  function  that  associates  a  number 
with  each  element  of  the  mesh,  which  represents  an  element 
contribution  to  the  global  error.  For  simplicity  we  choose 
an  error  indicator  based  on  the  jump  in  the  flux  —au'  of  a 


Fig.  8.  Cathode  reaction  rates  with  oxygen  limited  kinetics. 


scalar  u  between  adjacent  elements,  defined  on  an  element 
interface  v  by: 

Jx  =  —  a(x~)u\x~)  +  a(x+)u'(x+),  (22) 

where  the  minuses  (pluses)  denote  left  (right)  limits.  For  an 
element  K  =  (a,  b),  the  local  element  error  indicator  is  then 
defined  by: 

^-E^d^  +  l^bl2),  (23) 

where  the  summation  is  over  all  solution  components 
{(pQ,  (pm ,  X02}  and  over  all  data  points  (/;,  Vi). 

We  now  illustrate  the  effectiveness  on  the  estimation  of 
membrane  conductivity  amQ m  for  the  conventional  PEMFC 
data.  For  the  AMR,  we  sort  the  elements  by  their  error  indi¬ 
cator  values  and  choose  to  refine  those  elements  in  the  top 
20%  while  coarsening  those  elements  in  the  bottom  2%.  An 
element  is  refined  by  inserting  a  new  node  at  the  element 
midpoint,  while  a  pair  of  adjacent  elements  are  coarsened 
by  removing  the  common  node.  The  solution  on  the  original 
mesh  is  interpolated  onto  the  new  mesh  as  the  new  initial 
solution  guess  (for  implementation  see  [19]  or  for  general 
information  on  AMR  see  [20]).  In  Fig.  9  we  plot  the  error 
indicator  values  in  the  catalyst  layers,  where  the  error  indi¬ 
cators  take  on  the  largest  values.  We  can  clearly  see  that  the 
AMR  refines  the  mesh  in  order  to  equally  distribute  the  local 
error  indicator  values,  resulting  in  meshes  with  element  sizes 
varying  over  several  orders  of  magnitude.  The  result  is  that 
even  with  a  simple  error  indicator  we  can  reduce  the  num¬ 
ber  of  degrees  of  freedom  necessary  to  produce  a  reasonable 
parameter  estimate,  resulting  in  a  more  efficient  computation. 
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X  [1  e-4  m] 


x  [1  e-4  m]  Degrees  of  freedom 


Fig.  9.  Local  error  indicators  used  for  adaptive  mesh  refinement  in  the  anode  (top  left)  and  cathode  (top  right)  catalyst  layers.  Final  mesh  size  distribution 
(bottom  left)  and  convergence  of  membrane  conductivity  (bottom  right). 

5.  Conclusions  References 


In  this  study,  we  have  presented  an  algorithm  for  estima¬ 
tion  of  PEMFC  model  parameters  using  a  constrained  nonlin¬ 
ear  least  squares  algorithm.  Estimation  of  five  model  param¬ 
eters  in  a  simple  one-dimensional  electrochemistry  model 
using  two  different  experimental  polarization  curves  has  been 
demonstrated  using  this  procedure.  The  potential  for  further 
enhancing  the  algorithm  by  integrating  it  with  mesh  refine¬ 
ment  was  demonstrated  using  an  example  showing  the  greater 
efficiency  and  reliability  obtained  from  using  initial  guesses 
computed  on  coarse  grids,  as  well  as  locally  refined  meshes. 

This  work  naturally  suggests  potential  additional  research, 
including  attempting  parameter  estimation  for  more  complex 
PEMFC  model  equations,  two-  and  three-dimensional  sim¬ 
ulations,  as  well  as  transient  calculations.  A  key  point  is  the 
need  for  additional  data  besides  polarization  curves,  such  as 
local  data  along  the  channel  of  operating  cells,  net  water  trans¬ 
fer  balances  from  anode  to  cathode,  or  transient  polarization 
curve  data. 
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